home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 46 / Amiga Format CD46 (1999-10-20)(Future Publishing)(GB)[!][issue 1999-12].iso / -in_the_mag- / reader_requests / scilab / demos / arma / arma.dem
Text File  |  1999-09-16  |  2KB  |  84 lines
  1. mode(7)
  2. //
  3. // An example of arma simulation and identification 
  4. // form ( K.J. Astrom)
  5. // The armax process with the following characteristics 
  6. //    a=[1,-2.851,2.717,-0.865]
  7. //    b=[0,1,1,1]
  8. //    d=[1,0.7,0.2]
  9. // is simulated with an input u of a pseudo random binary type
  10. //
  11. // We use the simulated trajectory zd
  12. // as an input to the armax identification macro 
  13. // The noise in the armax is coloured so we can expect armax 
  14. // to give a biaised estimator 
  15.  
  16. a=[1,-2.851,2.717,-0.865]
  17. b=[0,1,1,1]
  18. d=[1,0.7,0.2]
  19. ar=armac(a,b,d,1,1,1);
  20. write(%io(2),"Simulation of an ARMAX process:");
  21. armap(ar);
  22. // The input
  23. u=-prbs_a(300,1,int([2.5,5,10,17.5,20,22,27,35]*100/12));
  24. // simulation with noise 
  25. zd=narsimul(a,b,d,1.0,u);
  26. // simulation without noise 
  27. z=narsimul(a,b,d,0.0,u);
  28. xselect();xbasc();
  29. plot2d2("enn",1,zd',[-1,1],"121","Simulated output");
  30. plot2d2("enn",1,1000*u',[-1,4],"100","Input [scaled]");
  31. halt();
  32. write(%io(2),"Identification ARX (least square):");
  33. [la,lb,sig,resid]=armax(3,3,zd,u,1,1);
  34.  
  35. // A bidimensional version of the preceding example 
  36.  
  37. a=[1,-2.851,2.717,-0.865].*.eye(2,2);
  38. b=[0,1,1,1].*.[1;1];
  39. d=[1,0.7,0.2].*.eye(2,2);
  40. sig=eye(2,2);
  41. ar=armac(a,b,d,2,1,sig);
  42. write(%io(2),"Simulation of an ARMAX process:");
  43. armap(ar);
  44. u=-prbs_a(300,1,int([2.5,5,10,17.5,20,22,27,35]*100/12));
  45. zd=narsimul(a,b,d,sig,u);
  46. z=narsimul(a,b,d,0.0*sig,u);
  47. write(%io(2),"Identification ARX (least square):");
  48.  
  49. [la,lb,sig,resid]=armax(3,3,zd,u,1,1);
  50.  
  51. // Spectral power estimation 
  52. // ( form Sawaragi et all) 
  53.  
  54. m=18
  55. a=[1,-1.3136,1.4401,-1.0919,+0.83527]
  56. b=[0.0,0.13137,0.023543,0.10775,0.03516]
  57. rand('normal');
  58. u=rand(1,1000);
  59. z=arsimul(a,b,[0],0,u);
  60. //----Using macro mese 
  61. [sm,fr]=mese(z,m);
  62. //----The theorical result 
  63. deff('[gx]=gxx(z)',['gx=(abs( a* exp(-%i*2*%pi*z*(0:4))'')**2)';...
  64.       'gx=abs( b* exp(-%i*2*%pi*z*(0:4))'')**2/gx']);
  65. res=[];
  66. for x=fr,res=[ res, gxx(x)];end;
  67. //----using armax estimation of order  (4,4)
  68. // it's a bit tricky because we are not supposed to know the order
  69. //
  70. [la,lb,sig,resid]=armax(4,4,z,u);
  71. a=la(1);
  72. b=lb(1);
  73. res1=[];
  74. for x=fr,res1=[ res1, gxx(x)];end;
  75. //
  76. leg="log(p) :using macro mese @ theoriqal value@log(p) : arma identifcation"
  77. xbasc();
  78. xtitle('Spectral power','frequency','spectral estimate')
  79. plot2d([fr;fr;fr]',[20*log(sm/sm(1))/log(10);...
  80.   20*log(res/res(1))/log(10);...
  81.   20*log(res1/res1(1))/log(10)]',...
  82.  [-2,-1,1],"111",leg, [0,-70,0.5,60]);
  83.  
  84.